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1.  Introduction 


In  practice  gravity  parameters,  such  as  geoid  heights,  or  deflections 
of  the  vertical  are  derived  either  by  numerical  integration  of  one  of  the 
well-known  integral  formulas  of  physical  geodesy,  or  by  an  optimal  estima¬ 
tion  technique  such  as  least-squares  collocation.  In  the  frame  of  the 
present  days  state-of-art--10"6  to  I0"7  relative  precision,  i.e.  e.g. 
geocentric  distances  to  about  0.5  to  5  meters— both  strategies  yield  reliable 
results  of  sufficient  accuracy.  But  in  view  of  both,  improving  measurement 
technology,  and  the  availability  of  a  rapidly  increasing  number  of  obser¬ 
vations,  one  aims,  mainly  in  the  context  of  geodynamics,  for  10"8  .  This 
high  goal  requires  however  not  only  improved  functional  models  but  at  the 
same  time  computation  techniques  that  make  adequate  use  of  the  physical 
information  contained  in  the  large  set  of  observations.  And  at  this  point 
the  individual  limitations  of  these  methods  become  visible. 

So  does  the  numerical  integration  approach  not  represent  an  estimation 
technique.  As  a  consequence  there  is  no  direct  way  to  associate  an  error 
measure  to  it.  Further  on,  not  the  original  (point)  observations  are  used 
for  the  integration  process,  but  block  mean  values  of  a  certain  block  size, 
derived  from  them.  But  there  exist  many  ways  to  compute  mean  values  and 
their  standard  deviations  from  the  point  observations.  In  addition,  one 
can  hardly  ever  be  sure  that  all  mean  values  entering  the  numerical  integra¬ 
tion  are  derived  by  the  same  processing  procedure. 

Least  squares  collocation,  on  the  other  hand,  does  work  with  the  original 
point  observations  and  is  an  optimal  estimator.  The  well-known  problem 
is,  however,  that  a  system  of  linear  equations  has  to  be  solved  of  a  size 
equal  to  the  number  of  observations.  Thus,  as  an  inherent  contradiction, 
the  more  data  are  collected,  or  the  smaller  the  sample  intervals  between 
them,  the  more  unstable  the  solution  becomes,  although  a  meaningful  global 
limit  (for  a  globally  continuous  data  coverage)  exists.  As  a  consequence, 
one  is  restricted  to  comparably  small  data  sets. 

Strategies  to  overcome  these  deficiencies  reach  from  (1)  the  arrange¬ 
ment  of  the  sample  points  in  a  regular  pattern  so  as  to  produce  certain 


favorable  matrix  patterns,  cf .  (Colombo,  1979),  via  (2)  sequential  estimation 
combining  potential  coefficients,  mean,  and  point  values,  cf.  (Tscherning, 
1975),  to  (3)  e.g.  a  combined  Integral  formula  and  collocation  approach, 
as  presented  in  (lachapelle,  (1977). 


The  method  presented  here  tries  to  combine  the  strong  points  of  the  nu¬ 
merical  integration  and  the  least-squares  collocation  method,  while  avoiding 
their  individual  deficiences.  It  is  based  on  the  global  limit  of  least- 
squares  collocation— where  an  analytical  inversion  of  the  in-the-global 
limit  infinite  dimensional  system  of  linear  equations  becomes  possible— 
and  maintains  therefore  its  character  as  an  optimal  estimator  with  an  error 
measure  associated  to  it.  This  global  limit  estimator  can  be  expressed 
as  a  stabilized  integral  formula.  A  numerical  integration  applied  to  it 
allows  to  process  arbitrarily  large  data  sets.  In  deriving  representative 
"area  weights"  it  becomes  possible  to  base  the  Integration  on  the  original 
point  values. 

The  two  essential  features  of  the  new  method,  the  derivation  of  a 
global,  stabilized  estimator  and  the  numerical  Integration  based  upon  the 
original  point  observations,  which  complement  each  other  in  an  ideal  manner, 
are  in  principle  independent.  So  could  one  apply  any  other  numerical  inte¬ 
gration  to  the  stabilized  integral  formula,  too,  e.g.  the  classical  one 
using  mean  block  values.  On  the  other  hand,  there  are  many  more  applications 
of  the  numerical  integration  based  on  point  values,  than  the  one  discussed 
in  this  context. 


2.  Stabilized  Global  Estimator 

In  (Neyman,  1974)  and  along  a  different  line  in  (Moritz,  1975)  a  so- 
called  "modified  Stokes  function"  Is  derived: 

sKtypij)  *  *  J2  nr  ^  pt(cos*pq)  (l) 

It  differs  from  the  Stokes  kernel  because  of  the  filter  factor 


In  equation  (1)  It  is 


Pj^cos^pg)... Legendre  polynomial  of  degree  l  , 
ij^pq  ...spherical  distance  between  P  and  Q  , 

cf  ...gravity  anomaly  model  degree  variance,  and 

<L  ...a  priori  error  degree  variance. 


We  shall  derive  in  the  sequel  this  equation  in  detail  for  an  arbitrary 
gravity  estimation  problem.  Let  us  assume  a  set  of  n  observed  ("~") 
gravity  parameters,  g-j  a  given.  As  unknown  the  global  disturbing  potential 
function,  t  ,  is  chosen.  Once  an  approximation  of  t  is  derived,  any  other 
functional  of  it  can  be  computed  straightforwardly.  The  character  of  t  as 
a  function  shall  be  indicated  by  expressing  it  as  a  vector  t  of  dimension 
(<»x  1).  To  be  more  precise,  we  assume  t  e  H  ,  a  Hilbert  space  with  reproducing 
kernel  K(P,Q) =  Ctt(P,Q)  .  The  spectral  components  or  coefficients,  t^  , 
of  t  form  an  infinite  but  countable  sequence,  or  vector  of  coefficients, 
eJt2  isomorphic  to  H  ,  see  e.g.  (Meissl,  1976).  It  is 


+i 


up)  *  i  i  if  t*1 1  ?..<p) . 

1=2  m=-l  rP  m  m 


(3) 


where  is  a  short-hand  notation  for  the  fully-normalized  surface  spherical 


harmonics  ofdegree  l  and  order  m 

P^mCsin^p)  cosmXp  for  ma 
P£ | m j  (sin4>p)  sin | m|Xp  for  m<  l 


Wp> 


(4) 


R  is  the  radius  of  a  convergence  (B jerhammar )  sphere,  and  rp  the  geocentric 
radius  of  P  .  For  K(P,Q)  we  choose,  cf.  (Moritz,  1980)  the  convergent 
series 


oo  2 

K(P,Q)  *  Ctt(P,Q)  -  l  (=V)a+1  c  P  (cos^p0)  ,  (5) 

with  the  disturbing  potential  degree  variances.  The  functional  relation¬ 
ship  between  the  observable  g  and  t  may  be  described  by  the  operator  equation 


g  ■  S_1  t  . 
nxl  nx»  »xl 
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(6) 


The  corresponding  relation  between  the  spectral  components  may  be  written  as 


9«m  =  Xl  4£m  • 

with  only  dependent  on  l  in  case  the  operator  S_1  is  isotropic,  cf. 
(Neyman,  1974).  For  a  detailed  derivation  in  terms  of  an  eigenvalue  decompo¬ 
sition,  see  (Rummel  et.  al.,  1979). 


Example:  In  case  g  is  a  vector  of  gravity  anomalies,  the  functional  is 

9  a  '  (f  +  3F  ) •  1  » 

the  so-called  fundamental  boundary  condition  of  physical  geodesy.  Together 
with  eq.  (3)  one  obtains  an  expression  of  the  type  of  eq.  (6) 


The  spectral  components  or  eigenvalues  of  S”1  are  in  this  case 


□ 

We  assume  E{§}  =  g  and  E{(g-g)(g-g)^ }  *  D.  The  least-squares  collocation 
estimate,  t  ,  ("A"  estimated)  is  derived  from  the  solution  of  the  variational 
problem 

min(JlS- +c*||  t  ||q-x)  ,  a>  1  ,  (8) 

cf.  (Rummel  et.  al.,  1979)  or  (Moritz,  1980,  ch.  28).  We  obtain. 

t  .  (S-V'S-  -  oC-)->  S’V1} 

«xl 

*  CttS'T  vS-xCttS-T  +  oD)'1  g 
®x°°  »xn  nx®  »x®  ®xn  nxn  nxl 


(9) 


The  terms  CS-T  and  S-1CS“'  in  (9)  express  the  application  of  the  operator 
S"1  to  ,  the  well-known  "propagation  law  of  covariances".  We  write 
in  short 

Cto  =  Qtt  S’T  •  and  (10) 

°°xn9 

C„  =  s-1c_s-T  .  (11) 


We  now  consider  the  global  limit  of  equation  (9)  where  §  becomes  a  random 
function  covering  the  entire  surface  of  the  earth  or  a  sphere,  e.g.  in  satellite 
altitude.  A  detailed  derivation  is  given  in  the  Appendix.  If  we  denote  the 
global  limit  of  the  linear  estimator  CttS-T(S-1CttS“T  +  aD)"1  by  lJ(P,Q) 
and  the  dimensionless  filter  or  smoothing  factor  in  eq.  (A-14)  by 

,  R,2(it+l)x2  = 

^  _  W  _  1 

ft  ■  ,r,2n+i),!=  . ■  t: •  U2) 


*0C  o  +  a°5  1  + 


'Alh2  ; 


We  obtain  from  equation  (A-14) 


‘-(VVW  *  ^  ^  (Q) 

-  I  MS?**1  ^  Mcos*pn)  • 


I  *'rp' 

For  our  example,  with  g  gravity  anomalies,  expression  (13)  becomes  the  mod¬ 
ified  Stokes  function  of  equation  (1).  It  would  be  identical  to  the  classical 
Stokes  function  for  fz-  1.  The  global  limit  of  the  optimal  estimator,  equation 
(9),  may  now  be  written  as 

t(P)  *  lim  {Ct _(C  +  aD)"1  g} 

dense  y 

-  lt  g 

«X“  “Xl 


Vi  I  tl  ««>  % 


a  % 


=  Tv  /  L(rp»rQ»^pq)  3W)  doq 


In  (Rummel ,  1982)  it  is  shown  that  equation  (14)  may  also  be  interpreted  as 
an  optimal  solution  of  the  stochastic  variant  of  the  global  geodetic  boundary 
value  problem,  with  only  one  realization  of  g  given  on  the  boundary  surface. 

The  product  X|  eg,  is  the  propagation  of  the  disturbing  potential  degree 
variances  to  the  degree  variance  of  the  observed  gravity  quantity.  As  to 
be  seen  from  equation  (12)  and  as  discussed  in  (Gerstl  &  Rummel,  1981)  the 
filter  factors  f^  depend  on 

-the  a  priori  error  variance  Og  , 

-the  regularization  factor  a  (usually  chosen  equal  to  1), 

-the  choice  of  the  radius,  R  of  the  Bjerhammar  sphere,  and 

-the  chosen  degree  variance  model  of  t  . 

The  stabilized  integral  kernel  L(rp,rg,^pg)  can  be  a  “tructed  very  easily 
for  any  problem  at  hand,  as  long  as  g  is  a  linear  fi  :ion  of  t  and  exists 
even  for  cases  where  the  non-stochastic  counterpart  do  ^ot.  It  is  e.g. 

possible  to  construct  a  global  "downward  continuation"  timator.  Since, 

because  of  ao,  ,  no  closed  formula  can  be  found  for  i  jn  (13),  it  is 

in  practice  approximated  by  the  Legendre  polynomial  truncated  at  a  very  high 

degree  Amax  . 

The  next  step  shall  be  to  prepare  for  a  finite  approximation  of  the  stabil 
ized  global  estimator  from  the  discrete  point  observations.  If  one  prefers 
to  work  instead  with  mean  values  of  a  certain  block  size,  the  traditional 
numerical  integration  techniques  may  be  directly  applied  to  equation  (14). 

3.  Area  Weights  for  Points  Irregularly  Distributed  on  a  Sphere 

The  difficulty  of  performing  a  numerical  integration  directly  based  on 
point  observations  irregularly  distributed  on  a  sphere  stems  from  the  problem 
that  area  weights  have  to  be  assigned  to  the  observations,  that  reflect  the 
local  data  density.  Area  weight  is  thereby  defined  as  a  unique  and  represen¬ 
tative  surface  area  element  (on  the  earth's  surface,  on  a  sphere,  e.g.  in 
satellite  altitude,  or  simply  on  a  unit  sphere)  assigned  to  an  observed  or 
derived  gravity  quantity  value.  When  working  with  mean  block  values,  they 
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are  e.g.  cos^AX  .  Since  the  area  weights  have  to  reflect  the  local  data 
density,  i.e.  small  area  weights  for  denser  point  samples,  one  has  first  to 
establish  a  quantitative  relation  between  the  points.  One  way  to  achieve 
this,  is  to  connect  the  points  by  a  triangular  net  with  minimal  side  lengths 
and  no  side  crossing  another  one,  as  displayed  in  Figure  1.  Then  one  could 
compute  the  area  of  each  triangle,  e.g.  MPQ,  and  associate  to  it  the  average 
of  the  observations  g.,-  at  the  three  nodes.  We  shall  follow  another  line: 

The  areas  of  all  triangles  adjacent  to  one  point,  e.q.  P  ,  compare  again 
Figure  1,  are  summed  up  and  divided  by  three,  because  each  triangle  is  shared 
by  three  points.  It  is  the  area  weight  wp  ,  assigned  to  the  discrete  observa 
tion  gp  at  P  .  The  weights  derived  this  way  are  unique,  reflect  the  local 
data  density,  and  their  sum  represents  exactly  the  total  area  covered  by  the 
measurements.  The  result  is  a  step  function  on  the  sphere  (or  the  surface 
of  the  earth)  with  step  area  equal  to  wp  and  step  size  (height)  equal  to 


One  could  think  of  many  sophistications  of  this  principle.  Instead  of 
this  discontinuous  step  configuration,  continuous  or  continuously  differential 
representations  could  be  designed  involving  more  and  more  of  the  neighboring 
observations  for  their  definition.  A  consequence  is  a  more  complicated  integra 
tion  procedure. 

The  problem  is  now  to  find,  first,  the  unique,  minimum  side  lengths  trian¬ 
gular  net  and,  second,  a  fast  algorithm  for  the  computation  of  the  area  weights 


3.1  Closest  Point  Triangulation  in  a  Plane 

The  problem  of  finding  a  unique,  non-overlapping  triangulation  to  irreg¬ 
ularly  distributed  points  ,  is  well-covered  in  the  literature,  cf.  (McCullagh 
&  Ross,  1980)  or  (Peucker  et.  al.,  1978).  We  shall  base  our  computations 
on  a  computer  algorithm  developed  by  M.  Gerstl ,  see  (Gerstl  et.  al . ,  1979), 
that  works  in  a  plane.  Its  theoretical  background  is  given  in  (Shamos  &  Hoey, 
1975).  The  concept  is  as  follows:  Given  is  a  list  of  Cartesian  coordinates 
and  yi  of  the  points  Qj  .  In  a  first  step  they  have  to  be  sorted  such 
that  xi  <  xi+i  and  y^  <  y-j+j  ,  if  xi  s  x^+j  .  Standard  procedures  per¬ 
forming  the  sorting  are  available  at  any  computer  installation.  The  first 
three  points,  Qx  ,  Q2  ,  and  Q3  of  the  sorted  list  form  a  triangle  and 
the  first  stage  of  a  convex  hull  that  shall  step  by  step  be  extended.  (K  is 
connected  with  those  vertices  of  the  convex  hull,  it  sees  (two  or  three). 

Any  side  facing  Q4  has  to  be  dropped  if  it  can  be  intersected  by  a  shorter 
one  originating  from  Q4  .  The  up-dated  convex  hull  contains  Q4  .  The  same 
principle  is  now  applied  to  every  subsequent  point  of  the  sorted  list.  This 
way  a  triangulation  is  spread  over  all  points,  with  no  intersecting  sides 
(except  at  the  nodes  Qi)  and  with  minimum  side  lengths.  As  a  final  product 
one  obtains,  first,  a  list  counter-clockwise  ordered  of  all  points  connected 
to  each  point  and,  second,  an  ordered  list  of  all  points  of  the  final  convex 
hull . 


Modifications  for  the  Work  on  a  Sphere:  The  same  principle  could  be 
directly  employed  to  data  given  on  a  sphere,  too.  But  in  order  to  maintain 
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the  efficiency  of  the  fast  computer  algorithm  in  this  case  one  has  to  avoid 
the  time-consuming  trigonometric  manipulations  on  the  sphere  as  much  as  possible. 
Therefore  we  compute  beforehand  from  the  given  latitude,  longitude,  and  radial 
distance  (or  constant,  radius)  3-dimensional  Cartesian  coordinates  of  all  points. 
Side-length  comparisons  can  then  be  based  on  chords  instead  of  spherical  dis¬ 
tances.  Azimuth  comparisons,  necessary  for  the  decision  which  sides  on  the 
convex  hull  are  facing  a  new  point,  are  carried  out  using  <&  and  X  as  if 
they  were  x  and  y  because  this  approximation  does  not  influence  the  compar- 
i  son . 


3.2  Data  Management  Aspect 

In  principle  the  triangulation  algorithm  could  be  applied  to  data  sets 
of  any  size.  But  the  program  has  to  use  an  integer  field  LC(LMAX.NMAX)  which 
contains  the  ordered  list  of  nodes  joint  to  each  point,  with  LMAX  the  total 
number  of  points  and  NMAX  the  maximum  number  of  joint  points.  Since  the  number 
of  joint  points  can  run  up  to  50  very  easily  at  an  intermediate  step  of  the 
processing  of  a  large  point  set,  one  arrives  with  point  sets  larger  than 
10,000  points  easily  at  the  storage  limits  of  even  very  large  computers.  The 
field  dimension  has  to  be  kept  so  large  because  it  is  at  least  in  theory  pos¬ 
sible  that  during  the  processing  changes  may  occur  even  in  the  innermost  zone 
of  the  convex  hull.  We  considered  up  to  now  three  possibilities  to  circumvent 
this  problem: 

1.  One  can  store  this  field  on  external  storage  and  work  with  direct  access. 

Then  the  computer  program,  which  takes  in  its  present  form  almost  no  c.p.u.- 
time  would  loose  much  of  its  efficiency. 

2.  Another  possibility  is  to  process  overlapping  point  sets  independently  and 
afterwards  merge  them.  Although  this  is  not  a  very  elegant  method,  it  works 
rather  efficiently  and  has  been  applied  for  the  numerical  example  following 
below. 

3.  One  could  work  with  the  maximum  available  storage  and  first  build  up  a  convex 
hull  in  this  limit.  During  the  processing  of  subsequent  points  there  exists 
theoretically  a  chance  that  sides  are  changed  in  the  inner  zone  of  the  convex 


hull.  But  the  chances  are  rather  slim.  Thus,  one  could  simply  freeze  the  Inner 
zone,  i.e.  not  allow  any  changes,  without  hardly  any  loss  in  the  criterion  of 
minimum  side  length.  The  frozen  part  could  be  stored  externally  and  the  cor¬ 
responding  core  storage  space  recycled  for  further  use.  This  possibility  is 
Investigated  at  the  moment. 


3.3  Area  Weight  Computation 


The  area,  F  ,  of  the  sum  of  m  spherical  triangles  adjacent  to  a  point 
P  ,  compare  Figure  2,  is: 
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with  F^  area  of  the  i-th  spherical  triangle  and  *aj+B^  +Yi  -180°  , 
the  corresponding  spherical  excess.  It  is  then 
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The  angles  are  obtained  from  the  spherical  azimuths  of  the  outside 
edges.  It  is 


*1  = 


”^Ai  ,i+l  “  Ai  ,i-l^  if  X1+l>X1>Xi-l 


(Ai,i+l  -  Ai,i-1)  otherwise. 


From  tanAi5i+!  and  tanAij.j  one  computes  tar<i  by  standard  trigonometric 
formulas.  In  order  to  avoid  the  evaluation  of  "arctan",  one  can  use  the  recursion 
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tan(  Z  <  ) 
i=l  1 


m-I 

tan(i|1  *.)  +  tan 
1-  (tan  mf1  <■)  tan< 


(17) 


for  the  computation  of  ^  .  Nevertheless,  all  spherical  azimuths  would 

have  to  be  computed. 

Again  an  approximation  is  Introduced  to  avoid  time-consuming  spherical 
computations.  Since  dense  data  coverage  is  the  underlying  motivation  for  the 
whole  method  the  type  of  approximation  used  here  is  justified.  We  use  the  Elling 
formula,  well-known  from  land  surveying,  and  find,  compare  Figure  3: 

Fp  1  i  cos*p  <♦,-♦,*,)(*,  +  Xm)(ji-)2  ,  (18) 

where  cos<j>p  takes  care  of  the  convergence  of  the  meridians.  For  a  block  en¬ 
closed  by  two  latitude  circles  and  two  meridians,  eq.  (18)  yields  the  usual 
area  element  cos<j>pA$AA.  No  trigonometric  function  have  to  be  computed  when 
using  eq.  (18).  Since  each  triangle  is  shared  by  three  points,  we  finally  find 
for  the  area  weight 

wp  "  tfp  *  -5-cos*p  X<V*i+i)(Vxi+i,(I55*)2  * 


(19) 


3.4  Final  Data  Base 


After  having  carried  out  these  computations  the  data  base  contains  besides 
the  original  entities,  e.g.  point  number,  coordinates,  measurement  value,  the 
computed  area  weight  wp  .  In  addition,  one  can  add  a  pointer  to  a  separate 
data  file,  that  contains  a  list  of  the  adjacent  nodes  to  each  point,  if  needed 
for  other  applications.  The  points  on  the  convex  hull  and  their  (distorted) 
area  weights  are  not  elements  of  the  data  base.  They  as  well  as  their  adjacent 
points  should  be  stored  in  a  separate  data  file.  This  way  it  is  possible  to 
restart  at  a  later  time  the  triangulation  of  a  new  area  leaving  the  old  data 
base  unchanged. 

New  discrete  points  Inside  the  original  point  set,  e.g.  new  measurements 
§  ,  can  be  Incorporated  In  such  a  way  that  only  that  triangle  Is  affected  In 
which  the  new  point  is  contained,  (compare  Figure  4).  This  way  only  four  area 
weights,  that  of  the  Incoming  point  and  those  of  the  three  vertices  of  the  tri¬ 
angle,  have  to  be  up-dated,  on  the  expense  that  locally  the  minimum  side  length 
criterion  is  probably  violated. 


Figure  4 


The  here  proposed  philosophy  of  working  with  the  point  observations  which 
are  related  through  a  "triangulation"  with  the  points  in  their  proximity  is 
supported  by  a  general  trend.  The  availability  of  very  powerful  data  manage¬ 
ment  systems  stimulates  to  work  with  the  large  sets  of  original  observations 
and  could  make  mean  block  averages  more  and  more  superfluous. 

Many  applications  in  geodesy,  geophysics,  surveying  and  civil  engineering 
can  be  envisioned,  such  as 

-  the  computation  of  gravity  parameters,  such  as  height  anomalies,  deflections 
of  the  vertical  etc.  from  large  sets  of  observations, 

-  the  estimation  of  surface  gravity  parameters  from  satellite  gravity  sensors 
(gradiometry  or  satellite-to-satellite  tracking), 

-  topographic  or  topographic-isostatic  reduction  and  inverse  geophysical  modell¬ 
ing  from  digital  terrain  models,  or 

-  height  interpolation,  mass-volume  estimation,  profile  computation  etc.  from 
photogrammetric  or  tachymetric  digital  terrain  models. 


4.  Error  Estimate 


As  important  as  the  estimate  t  itself  is,  we  must  have  a  reliable  error 
measure  for  it.  For  a  moment,  we  still  remain  with  the  global  limit  case,  where 
§  is  assumed  to  be  a  random  function.  The  posteriori  variance-covariance  matrix 
of  t  eq.  (9),  resulting  from  least-squares  collocation  is,  cf.  (Moritz, 1980) : 

-tt  *  -tt  ‘  -tt-  T  (-”1 -tt-  T  +ct3)"1  S'1  ?t1 
aovao  ooxoo  «x ooooxn  ®nx«  oo  x»  ooxn  nxn  nx® 00 x00 


(20) 


(21) 


or  with  the  defined  linear  estimator,  eq.  (13), 

i«  ■  •  tTr  ctt . 

®x»  ®x»  “xn  nx»  ®x® 

From  a  simple  manipulation  we  obtain 

§tt  =  Qtt  •  -T  {(§'x  5tt  S"T  +  “Q)  (§"x  9tt  §"T  +  “S)’1}  # 


a  Ctt  -  LT  (p  Ctt  S‘T  +  oO)  L  (22) 

Furtheron  with  a  similar  manipulation  one  derives  from  equation  (201: 


-tt  =  9tt  "  9tt  §  T  L  -  LT  §-1Ctt  +  LT(S*lCtt  S'T  +  oO)  L 
-  (LTS~1  -  I)Ctt  (Js-1-!)7  +  aLTD L 

*  BTCttB  +  aLT  0  L  (23) 
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with  the  unit  matrix  I  and  the  bias  expression 

Bt  =  I  -  LT  S"1 
«x«  ®x®  ®xn  nx® 


(24) 


Equations  (21),  (22),  and  (23)  are  three  alternative  forms  to  equation  (20). 
Equation  (22)  Is  convenient  from  the  numerical  evaluation  point  of  view,  because 
the  inversion  of  the  full  matrix  (§-1CttS"T+ aO)  is  avoided,  once  L  is  known. 
Equation  (23)  is  Interesting  from  the  interpretation  point  of  view.  Equations 
(20)  to  (24)  are  still  infinite  dimensional.  For  the  evaluation  of  the  a  poster 
iori  variance  mp  at  P  ,  (or  analogously,  the  covariance  between  two  arbitrary 
points  P  and  Q)  it  is  convenient  to  introduce  the  evaluation  functional  ep  . 
We  define  e.g. 
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Then  equations  (22)  and  (23)  become 


4  ■  Ctt,P  -  ij  (§-lCt1.S-T  ♦  CD)  Lp 

‘  Ctt,P  '  tp  <£gg  *  “9)  tp 

(26) 

and 

"•p  *  §P  5tt  -P  +  a-P  5  kp  • 

(27) 

With  equations  (A-3)  and  (A-14)  of  the  appendix  we  find  for  the  spectral  form 

of  equation  (26): 
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where  we  used  that  £  ^m(P)  Y*m(P)  *  (2£+l)  P^l)  »  2t+l  . 

From  equation  (12) 

one  obtains  for 
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(29) 

Together  with  eq.  (28)  it  yields  the  limits 
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l.e.the  a  posteri  variance  becomes  equal  to  the  a  priori  one.  All  these  error 
estimates  hold  for  the  case  where  a  global  random  function  §  is  assumed  to 
be  given.  In  reality  we  assume  densely  spaced  point  observations  ,  i«l,...n  . 

Equation  (14)  is  approximated  by  a  numerical  integration  procedure.  Equation  (14) 
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is  replaced  with  very  little  loss  by 

*p  *  i^1  L ’ ^rp»rQ|^PQ1 ^  h  ttV  •  ^ 

In  equation  (30)  the  w.j  are  the  area  weights  as  obtained  from  eq.  (19),  and 
L‘  is  the  discrete  approximation  of  the  operator  L  ,  eq.  (13).  Important 
is  thereby  that  for  the  computation  of  L*  ,  a§  has  to  be  replaced  by 

d*  »  al  (A- 10) 

t  n 

with  as  ■  ~  l  wi  ,  the  average  area  weight.  Thus,  the  a  posteriori  variance, 
equation  (27)  jBicomes 

m£  *  Cu  Bp  ♦  aLpT  0  Lp  (31) 

def 

i.e.  in  all  expressions  L  is  replaced  by  L'  (a  is  usually  chosen  a  *  1). 

The  second  term  on  the  right-hand  side  of  equation  (31)  becomes  for  uncorrelated 
observations  with  homogeneous  a  priori  variance 

°o  tp...  the  pure  error  propagation.  This  is  the  only  part  usually 
considered  in  geodetic  numerical  Integration  procedures.  It  is  also  well-known 
from  the  least-squares  adjustment  of  overdetermined  problems.  As  Important 
is  the  first  term 


Bp  Ctt  Bp...  the  discretisation  error.  From  the  definition  of  B  , 
equation  (24),  one  sees  that  this  part  becomes  zero  for  global  and  continuous 
data,  free  of  errors.  For  a  large  amount  of  observations.  Its  evaluation  is 
rather  tedious,  because  of  the  size  of  Ct4.  . 


5.  Numerical  Example 

In  a  first  and  still  preliminary  test  the  method  was  tested  through  a  com¬ 
putation  of  some  l°xl°mean  gravity  anomalies  from  adjusted  GEOS-3  altimeter 
data  in  the  North  Atlantic.  This  example  was  chosen  for  several  reasons: 

-  The  GEOS-3  altimeter  data  fulfills  very  well  the  requirement  of  being  at  least 
locally  dense  and  homogeneous. 

-  The  computation  of  gravity  anomalies  from  the  "quasi"  geoid  undulations  as 
obtained  from  altimetry,  poses  a  very  unstable  problem.  Thus,  it  is  a  good 
test  for  the  stability  behavior  of  this  numerical  integration  approach. 

-  In  the  same  geographical  area  l°x  1°  anomalies  are  available  as  obtained  from 
shipborne  gravity  measurements,  as  well  as  those  estimated  from  GEOS-3  altimetry 
by  least-squares  collocation,  compare  (Rapp,  1977).  With  collocation  typically 
around  400  observations  could  be  processed  per  gravity  anomaly  estimate  out 

of  the  total  number  of  500,000,  whereas  with  the  new  method  we  could  easily 
work  with  7000  without  even  taking  into  account  the  data  management  aspects 
mentioned  in  Chapter  3.2. 


In  this  test  computation  approximately  7000  adjusted  sea  surface  heights 
as  derived  from  GEOS-3  altimetry  were  included.  At  the  accuracy  level  of 
GEOS-3  (±0.5  to  0.8  m)  they  can  be  considered  synonymous  to  geoid  undulations. 

The  adjustment  of  the  altimeter  data  is  described  in  (ibid.).  The  processing 
steps  of  the  new  method  are  displayed  in  a  flow  chart  in  Figure  5.  First, 
for  all  points  the  area  weights  were  computed  according  to  eq.  (19)  and  included 
in  the  data  file.  Prerequirement  was  the  determination  of  a  triangulation 
connecting  all  points,  as  described  in  Chapter  3.1  The  triangle  net  is  shown 
in  Figure  6,  where  one  can  still  recognize  the  original  pattern  of  the  satellite 
ground  tracks.  In  an  enlarged  sub-area  the  principle  of  the  area  weight 
computation  is  illustrated,  Figure  7.  Independent  thereof,  a  table  with 
the  function  values  of  the  stabilitzed  Integral  kernel,  equation  (13),  was  computed 
at  small  Intervals  A*  -  1\  i.e.  at  (*k  |0°  (1*)  15°)  .  During  the  numerical 
integration  process  it  is  linearily  interpolated  between  the  tabulated  values. 

The  computation  of  the  table  values,  theoretically  a  summation  up  to  infinity, 
has  been  truncated  at  £max  -  10,000.  The  degree  variance  model  by  Tscherning 
and  Rapp  (1974)  was  chosen  for  the  definition  of  the  c‘t  ,  compare  again  eq.  (13). 
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For  the  a  priori  error  variance,  o§  ,  we  chose  0.64  m2 (std.  dev.  ±0.8  m). 
According  to  equation  (A-10)  o§  was  multiplied  by  ,  with  As  the  average 
area  weight.  The  obtained  d'  is  used  instead  of  ao§  in  eq.  (12).  The  tab¬ 
ulated  kernel  values  together  with  the  area  weights  and  the  altimeter  derived 
geoid  undulations  allow  to  perform  the  numerical  integration,  equation  (30), 
completely  along  the  line  of  e.g.  the  numerical  Integration  of  the  Stokes'  inte¬ 
gral  In  the  context  of  geoid  computations.  The  long-wavelength  part  of  the 
undulations  and  of  the  mean  gravity  anomalies  to  be  estimated  has  been  treated 
separately  by  subtracting  from  the  altimeter  derived  undulations  the  contri¬ 
butions  coming  from  a  set  of  potential  coefficients  and  adding  back  the  cor¬ 
responding  contribution  to  the  mean  anomalies.  For  the  long-wavelength  part 
the  Goddard  Earth  Model  9  set  of  potential  coefficients  up  to  degree  20  has  been 
used.  Consequently  the  summation  of  the  stabilized  integral  kernel,  eq.  (13), 
has  been  started  at  degree  21.  Table  1  shows  the  results  for  a  set  of  arbitarily 
selected  l°xl°  blocks.  Column  three  gives  the  mean  anomaly  values  as  obtained 
from  the  GEM-9  model,  column  four  the  terrestrial  (shipborne  gravimetry)  mean 
anomalies,  column  five  the  values  obtained  with  the  new  method,  column  six  the 
collocation  results  provided  by  Or.  Rapp.  The  estimates  derived  with  the  new 
approach  show  a  good  agreement  with  the  terrestrial  and  collocation  data  with 
an  overall  tendency  to  be  closer  to  the  terrestrial  anomalies.  It  would  be 
unreasonable  to  draw  any  further  reaching  conclusions  from  this  limited  test. 

This  application  can  be  considered  an  alternative  to  the  techniques  described 
in  (Rummel  et.  al.,  1977). 
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Table  1;  A  comparison  of  rxl°mean  gravity  anomalies  estimated  from  geoid 
undulations  (derived  from  GE0S-3  altimetry)  by  least  squares  collocation  (column 
6)  and  by  the  new  technique  (5)  with  the  corresponding  terrestrial  values  (4). 
The  mean  anomalies  as  obtained  from  the  GEM-9  set  of  potential  coefficients 
up  to  s,-=20  is  contained  in  column  3.  The  coordinates  $  and  A  of  the  north¬ 
west  corner  of  each  block  are  given  in  columns  1  and  2.  (Dimension  Ag:  mgal  * 
10"5  ms'2) 


6.  Further  Sophistications 

One  of  the  most  important  features  of  least-squares  collocation  is,  that 
it  allows  to  combine  in  an  optimal  way  heterogeneous  data,  i.e.  different  gravity 
quantities  distributed  In  an  arbitrary  manner  over  various  points.  Since  the 
present  method  is  derived  from  a  global  limit  process,  this  feature  is  lost 
in  this  generality.  With  other  words,  there  is  no  simple  way  to  combine  a  gravity 
anomaly  at  one  point  with  a  deflection  of  the  vertical  component  at  another 
one,  and  a  potential  difference  at  a  third  point.  Sjbberg  (1979)  has  shown, 
however,  that  It  Is  possible  to  combine  heterogeneous  data  in  the  global  limit, 
i.e.  it  is  again  possible  to  derive  a  global  limit  estimator  analogous  to  chapter 


data  file 


numerical  integration 
"gravity  parameter  estimation 


2,  that  combines  global  random  functions  of  different  type.  In  its  most  gen¬ 
eral  form  the  optimal  estimator  may  be  written,  analogous  to  eq.  (14),  as: 


with  m  different  gravity  quantities  g^J^  ,  j=l,...m.  The  expression  for 
the  vector  estimator  becomes,  analogous  to  eq.  (13): 

-4rP’rQ’W  (jFj^*+1(2l+l)  P*(cos*pq)  (33) 
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In  eq.  (33)  is  the  diagonal  matrix 
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where  the  X 
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servable  quantities  g 
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are  the  spectral  coefficients  relating  the  ob- 
(Q)  to  t(P),  compare  eq.  (7).  The  dimensionless 


filter  vector  becomes,  analogous  to  eq.  (12): 
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where  s  -  (- — )2  and  all  regularization  factors 
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were  chosen  equal 


to  one.  For^the  case  of  only  two  different  observed  gravity  quantiites,  g^ 
and  ai2^ 


the  inversion  contained  in  eq.  (35)  can  be  performed  explicitly. 


We  find  after  inversion  (A. . .determinant  of  the  matrix): 


Equations  (32)  to  (35)  are  of  practical  relevance  in  cases  where  the  same  com¬ 
bination  of  gravity  quantities  is  available  in  all  of  the  densely  spaced  but 
discrete  sampling  points.  Again  the  same  numerical  integration  procedure  can 
be  employed.  A  typical  example  would  be  a  satellite  gradiometer  mission,  where 
in  an  ideal  case  in  each  sampling  point  five  linearily  independent  second-order 
derivative  components  of  the  gravitational  potential  are  observable.  The  bene¬ 
fit  of  having  per  point  five  components  instead  of  only  one  can  easily  be  inves¬ 
tigated  with  the  formula  apparatus  presented  here. 

The  new  method  should  permit  a  better  exploitation  of  the  information  content 
of  large  and  densely  spaced  sets  of  observed  gravity  quantities,  prerequirement 
for  a  step  towards  a  10"8  relative  accuracy.  We  did  not  at  all  go  into  the 
question  what  model  improvements  are  required  at  the  same  time  to  achieve  this 
goal.  But  since  the  new  method  is  set  up  the  same  way  as  the  classical  numerical 
integration  techniques  in  physical  geodesy,  all  model  sophistications  developed 
for  them  can  be  applied  here,  too.  This  includes  especially  ellipsoidal  cor¬ 
rection,  and  tidal,  topographic  and  ellipsoidal  reduction,  as  discussed  in  the 
literature,  cf.  (Mather,  1973),  (Moritz,  1974),  (Rummel  &  Rapp,  1976),  (Rapp, 
1981),  (Groten,  1982),  although  on  the  other  hand  the  reduction  concept  is  not 
very  attractive  from  the  theoretical  point  of  view,  and  should  eventually  be 
replaced  by  more  sophisticated  models.  Also  the  numerical  side  can  be  improved, 
e.g.  by  minimizing  the  truncation  error,  when  working  with  data  limited  to  a 
certain  cap,  compare  (Jekeli,  1980),  or  by  optimizing  the  combination  with  sets 
of  potential  coefficients,  as  proposed  by  Wenzel  (1981). 


7.  Conclusions 


From  the  application  point  of  view,  the  classical  numerical  integration 
procedures  in  physical  geodesy  have  two  disadvantages.  They  are  not  estimators 
and  are  therefore  not  linked  to  a  stochastic  model  and  they  work  with  mean  block 
values,  which  come  from  some  subjectively  chosen  pre-processing.  Linear  estima¬ 
tion  techniques,  such  as  least-squares  collocation,  on  the  other  hand,  who  do 
not  have  these  drawbacks,  require  usually  the  solution  of  a  large  system  of 
linear  equations,  which  limits  their  applicability.  The  technique  presented 
here  tries  to  avoid  these  problems.  It  represents  a  "global  limit"  estimator. 
Especially  when  treating  unstable  problems  this  offers  considerable  advantages 
against  the  classical  integral  procedures.  In  addition,  the  error  measure  does 
not  contain  only  the  pure  error  propagation,  but  also  the  discretisation  error. 
Independent  of  these  features  the  numerical  integration  is  based  on  discrete  point 
values.  We  hope  to  initiate  a  trend  in  this  direction,  which  seems  to  us  timely, 
considering  the  fantastic  prospects  of  modern  data  management  capabilities.  The 
efficient  processing  of  even  very  large  sets  of  original  observations  should  be 
no  problem  in  the  near  future.  The  numerical  integration  method  itself,  and 
the  underlying  triangulation  of  "digital  terrain  model"  type  of  data  sets  offers 
many  more  applications  than  the  one  presented  in  our  example,  in  geodesy,  as  well 
as  in  surveying,  photogrammetry ,  civil  engineering,  and  geophysics.  Examples 
and  envisioned  improvements  to  be  implemented  in  the  numerical  procedure  are 
lined  out  in  Chapter  3. 

The  also  described  capability  of  optimal  data  combination  based  on  the 
ideas  of  L.  Sjdberg,  could  be  especially  useful  for  the  study  of  future  dedicated 
satellite  gravity  missions,  to  get  an  insight  into  the  improvements  coming  from 
simultaneously  sensing  the  gravity  field  in  different  spatial  directions. 

Necessary  model  improvements,  on  the  one  hand,  and  optimal  data  processing, 
on  the  other  hand,  seem  to  us  the  major  obstacles  on  the  way  to  a  " 10" 8 -precise 
gravimetric  geodesy".  This  study  is  supposed  to  be  a  contribution  to  the  latter. 
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Appendix 


The  stabilized  global  estimator  is  obtained  by  considering  the  limit  case 

i 

of  equation  (9),  where  the  observation  points  become  globally  dense,  i.e.  where 
one  approaches  a  continuous  coverage  with  observations  all  over  the  earth.  This 
type  of  limit  is  explained  in  (Moritz,  1975,  ch.2).  The  derivation  of  the  limits 
is  especially  convenient,  if  one  expresses  all  quantities  in  their  spectral 
form.  For  then,  as  pointed  out  in  chapter  2,  the  functions--or  global  limits- 
can  be  represented  as  infinite  dimensional  vectors.  The  inner  product 
(f  ,  g)H  in  H  for  example  defined  as 

If  .  g|„  *  jV  /  f(P)  g(p)  <iop 

a 

becomes,  when  using  the  spectral  forms. 


{f  •  9>h  -  /  O  hm  ’im  <p» 
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l  f£m  9tm  ’ 


tm 


where  the  orthogonality  relationships  of  spherical  harmonics  were  used.  In  a 
first  step,  we  shall  now  apply  this  principle  to  equations  (10)  and  (11),  where 
the  outcome  is  known  from  the  literature.  The  (isotropic)  operator  S"1  can 
be  written  as 


S'1 

®xn 


I  <7^  h  (cos*Q'Qjl 

1*1 


where  the  addition  theorem  was  used: 


.n. 


(A— 1 ) 


^(cosi^q.q)  =  2^  l  7£m(Q)  . 


(A-2). 


Applying  the  principle  shown  above  we  find  with  eq.  (4)  for  eq.  (10): 
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Furtheron,  for  equation  (11)  one  obtains: 
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(A-4) 


Now  to  the  global  limit,  which  we  denote  1im(Qj  globally  dense  on  Oq)  or 
short  g^ense  •  In  case  random,  uncorrelated  measurement  noise  with  variance 
ai  Is  assumed  the  global  limit  of  the  a  priori  variance-covariance  matrix  0 
can  be  expressed  as 


11m 

dense 


D-D 

O »X« 
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compare  (Helskanen  &  Moritz,  1967,  ch.  7-7).  It  Is  in  principle  a  critical, 
point,  because  several  severe  statistical  assumptions  enter  this  derivation. 

We  try  a  short  derivation.  Let  us  assume  the  a  priori  model  of  the  uncorrelated 
noise  Is  expressable  by  the  random  function  2( P)  .  Its  covariance  function  is 


D(P,Q)  -  ai  6 (P ,Q) 


(A-6) 


with  the  delta  function 

I«  for  P  »  Q 

and  P.Q  e  o 

0  for  P  f  Q 

The  spectral  components  of  D(P,Q)  (power  spectrum  components,  in  the  time 
series  terminology)  are  obtained  from: 


D(P.Q)  »W(PQ)  VQ)  d°P  d°'q 
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The  mathematical  distribution,  5(P,Q)  ,  Is  defined  in  such  a  way  that 

■fg  /«<P.P>  Yu  (?)  d<Jpd«f-l  (A-8) 

o 

Thus,  one  obtains  from  equation  (A-7) 


d«m  "  d  *  •  (A-9) 

For  the  finite  case  with  one  sample  point  per  As  =  sin  4>A4>AA  on  the  unit  sphere 
.o  ,  Jekell  and  Rapp  (1980)  derived  the  approximation  to  equation  (A-9): 

d'  *  -  o?  .  (A- 10) 

With  equations  (A-4)  and  (A-9)  the  global  limit  of  (S"lCttS“T-»aD) ,  compare  equation 
equation  (9),  becomes 

dense  §'T 
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where 

The  inverse  is 
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Finally,  we  obtain  for  the  global  limit  of  the  complete  estimator  with  equation 
(A-3): 
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